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Abstract—We introduce a new finite-difference time-domain 
(FDTD) method for solving the TM (transverse magnetic) mode 
reduction of Maxwell’s equations with complicated material 
interfaces. The method uses an unstructured quadrilateral mesh 
intended to conform to the complicated geometry of material 
interfaces. Using a quadrilateral to unit-square local coordinate 
transformation, we remap all the field components from the 
physical grid to a computational grid consisting of unit-squares. 
On the computational grid of unit-squares we update the fields 
in time using the usual FDTD field updates. The local coordinate 
remap, however, changes a scalar material on the physical space 
to an anisotropic material on computational space; we handle the 
mixed components of the anisotropic material by interpolating 
the non-colocated field components. Finally, we resolve material 
parameter jumps across interfaces using a harmonic-mean. By 
testing our method on the TMo,; mode in partially filled infinite 
cylindrical cavity, we verify that the method converges, is second 
order accurate, and maintains stability. 


I. INTRODUCTION 


The FDTD method [1] is one the most popular compu- 
tational methods used for solving electromagnetic problems. 
Rightfully so, it has many desirable properties such as explicit 
time-stepping and simple parallelization. However, FDTD 
does suffer from key problem dependent issues; the most 
pertinent of which is staircasing of boundary conditions and 
material interfaces. The main tools to address staircasing 
in FDTD are methods which conform the grid to perfect 
electric conductor (PEC) boundary conditions [2] and handle 
interfaces with sub-cell averaging techniques [3]. 

Nonetheless, these fixes skirt the crux of the issue that 
structured rectangular grids do not conform to complicated 
geometry. To address it directly, researchers have developed 
several unstructured grid approaches which use FDTD on 
conformal grids produced by meshing tools. These include un- 
structured triangle-polygonal meshes [4] and non-orthogonal 
mixed-polyhedral unstructured meshes [5]; all of which per- 
form the FDTD update on the physical grid (as opposed to a 
computational grid). 

Contrasting with prior unstructured conforming mesh ap- 
proaches, our approach borrows a technique from the finite- 
element community [6] in which a physical grid element is 
mapped to the unit-square where time updates will take place. 
Then we transform back to the physical grid to measure fields. 

For this conference paper we describe Maxwell’s equations 
under a coordinate transformation, introduce the unit-square to 
quadrilateral transformation (called the Piola transformation 
that’s equivalent to transformation optics), describe our TM 
mode equation FDTD method on an unstructured quadrilateral 
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Fig. 1. Depiction of the unit-square to quadrilateral bilinear transformation 
used in unstructured mesh FDTD method. Equivalent Faraday law integration 
paths for physical and computational cell are shown in red. 


grid, and verify our method. We finalize by summarizing the 
results in a conclusion. 


Il. MAXWELL’ S EQUATIONS UNDER A COORDINATE 
TRANSFORMATION 


Maxwell’s equations 


aD=VxH, 3B=-VxE (1) 


with the constitutive relations H = p~'B and E = e~'D 
are invariant under coordinate transformations. That is, under 
a map f : R? — R?® with Jacobian 7 = [0, f,0,f,0-f], the 
coordinate transformed fields E = 7~7E and H = J~7H 
satisfy Maxwell’s equations with the transformed constitutive 
parameters #2 = cu I and € ae a . Note that the 
tranformation may make 7 and € anisotropic. 

For the TM mode equations (0, = 0, z = Oand H, = FE, = 
E,, = 0) consisting only of field components (H,, H,, E.), 
the transformed fields and constitutive parameters become 

H=7 'H, E,=E. 

| nee 

~ det J’ iam det 7° 

where 7 is now a 2 x 2 matrix. Using this we may apply the 
inverse of the Piola transformation to remap every cell (each of 
which has its own local transformation) from an unstructured 
quadrilateral mesh to a computational grid consisting only of 
unit-squares. The bilinear transformation together with Piola 
transformation, depicted in Fig. 1, is given by 


(2) 


(3) 
(4) 
where M = P, — P; — (P3 — Po). Once on the computational 


grid of unit-squares, we apply an FDTD update, the details of 
which are described in the next section. 


f(x,y) = Po + (Pi — Po)z + (Ps — Po)y + mary 
J(w.y) = [Pi — Po + My, P3 — Pp + Mz], 
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Il]. TM MODE FDTD METHOD ON AN UNSTRUCTURED 
QUADRILATERAL GRID 


Our TM Mode FDTD method places FE, at quadrilateral 
midpoint (average of all corners) and places H tangential 
components at edge midpoints (continuity of Phisiicng)s see 
Fig. 1 for a depiction of field component locations. In order 
to consistently orient the edges of the quadrilateral grid (so 
that we may compute the circulation of H fields), we use the 
algorithm described in [7]. Then we map the fields and ma- 
terial parameters from the physical grid to the computational 
grid using the inverse of the transform in (2) with Jacobian 
(4) WH = JH, Ey = Ey op = Ceaty)s pr" and 
€ = €det /. On the unit-square grid we update the fields with 
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D,\*) = D,/|* + At) oH? 


(5) 

i=0 
Bee = Bee a3 At(E,|{** = EG") (6) 
Rae Dy, Hau. 8 (7) 


where o; is the orientation the square’s edge and the index | 
refers to the square on the positive side of the edge’s orien- 
tation; 0 means the negative side. The effective permeability 
is given by the harmonic average pj, = (uo + z+). In 
order to handle the anisotropic part of ihe we interpolate non- 
colocated fields to their required location. For example, to get 
Ay = Wy Be + Wy By we interpolate B, to the location of 
B,,. Once the time-updates are finished we switch back the to 
physical space using the forward transform in (2). 


IV. METHOD VERIFICATION 


We verify the method converges, is second order accurate, 
and is stable by testing it on a TMo, mode solution in a 
partially filled infinite cylindrical cavity with perfect magnetic 
conductor (PMC) boundaries. The cavity of interest with a 
sample quadrilateral mesh, generated by gmsh [8], is depicted 
in Fig. 2. An exact solution for FE, in the cavity is given by 


E%(p, 6, t) = [A*Jo(B’p) + B’Yo(8'p)|cos(wt) (8) 


where Jo and Yo are the O order Bessel functions of the Ist 


kind and 2nd kind, respectively. The constants A’, B’, 6°, 
Bi 


Ww = = are determined by ,4;, the interface condition, and 
the PMC boundary condition [9]. By varying the mesh size 
and comparing the numeric solution to the exact solution in 
the L? norm at time t = 3, we find that the method converges 
with 2nd order accuracy, results are plotted in Fig. 3 (a). 
Furthermore, we compute the eigenvalues of our method for 
the mesh in Fig. 2 and plot them in Fig. 3 (b). We find they 
all lie on the unit circle, confirming stability for the TMo1 


problem. 


V. CONCLUSION 


We introduced an unstructured mesh FDTD method that 
uses a local coordinate transformation to perform field updates 
on a computational grid consisting of unit-squares. We then 
verify it with an infinite cylindrical cavity that is filled with 
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Fig. 2. Mesh of infinite cylindrical cavity with e = 1 and PMC boundary 
at p = 1. For 0 < p < 1/2 (gray shaded area) it is filled zo = 1 material. 
For 1/2 < p <1, it is filled py = 2. 
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Fig. 3. Plot (a) shows convergence of the TMo,1 solution to partially filled 


cylindrical cavity. The observed L? error for several mesh sizes are plotted 
along with a reference 2nd order convergence line in blue. Plot (b) shows the 
method’s eigenvalues for the mesh in Fig. 2; all lie on the unit circle. 


two different permeabilities. We find that it converges, is 2nd 
order accurate, and is numerically stable. For future work we 
will investigate adding frequency dependent material to the 
method and extending it to three-dimensions. 
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